View surface for a subject or group, overlay brain mask, labels, and statistics.
In [1]:
"""Make static images of lyman results using PySurfer."""
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import os
import os.path as op
import sys
import tempfile
import argparse
from textwrap import dedent
from time import sleep
import numpy as np
from scipy import stats
import nibabel as nib
import matplotlib.image as mplimg
from mayavi import mlab
from surfer import Brain, io
import surfer as sf
import moss
import lyman
%gui qt
[-h] [-subjects [SUBJECTS [SUBJECTS ...]]] [-experiment EXPERIMENT] [-altmodel ALTMODEL] [-level {subject,group}] [-regspace {mni,fsaverage,epi}] [-output OUTPUT] [-geometry GEOMETRY] [-nowindow] [-keepopen]
In [113]:
class Args(object):
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
args = Args(subjects='~/Experiments/ObjFam/data/subids_subset_no23or19.txt',
experiment = 'objfam',
level = 'group',
regspace = 'fsaverage',
geometry = 'inflated',
output = 'group',
altmodel = 'studyrep_hpf256')
# altmodel = 'morphmem_schacter')
In [114]:
config_opts = dict(background="white", width=600, height=370)
In [115]:
# Load the lyman data
subjects = lyman.determine_subjects(args.subjects)
project = lyman.gather_project_info()
exp = lyman.gather_experiment_info(args.experiment, args.altmodel)
contrasts = exp["contrast_names"]
z_thresh = exp["cluster_zthresh"]
# Get the full correct name for the experiment
if args.experiment is None:
exp_name = project["default_exp"]
else:
exp_name = args.experiment
exp_base = exp_name
if args.altmodel is not None:
exp_name = "-".join([exp_base, args.altmodel])
In [116]:
# Group-level
# ===========
if args.level == "group":
temp_base = op.join(project["analysis_dir"], exp_name, args.output,
args.regspace, "{contrast}")
subj = 'fsaverage'
if args.regspace == "fsaverage":
sig_thresh = -np.log10(stats.norm.sf(z_thresh))
sig_thresh = np.round(sig_thresh) * 10
corr_sign = exp["surf_corr_sign"]
sig_name = "cache.th%d.%s.sig.masked.mgh" % (sig_thresh, corr_sign)
stat_temp = op.join(temp_base, "{hemi}/osgm", sig_name)
mask_temp = op.join(temp_base, "{hemi}/mask.mgh")
png_temp = op.join(temp_base, "{hemi}/osgm/zstat_threshold.png")
png_colorbar = op.join(temp_base, "{hemi}/osgm/colorbar.png")
png_label_temp = op.join(temp_base, "{hemi}/osgm/zstat_threshold_{label}.png")
png_roi_temp = op.join(temp_base, "{hemi}/osgm/zstat_threshold_{roi}.png")
else:
stat_temp = op.join(temp_base, "{hemi}.zstat1_threshold.mgz")
mask_temp = op.join(temp_base, "{hemi}.group_mask.mgz")
png_temp = op.join(temp_base, "zstat1_threshold_surf.png")
png_label_temp = op.join(temp_base, "zstat_threshold_{label}.png")
corr_sign = "pos"
In [117]:
stat_temp
Out[117]:
In [7]:
def plot_colorbar(zthresh, zmax, cmap, filepath):
fig = plt.figure(figsize=(10,.3))
cbar_height = 1
cbar_ax = fig.add_axes([.3, .01, .4, cbar_height - .01])
cbar_ax.set(xticks=[], yticks=[])
bar_data = np.linspace(0, 1, 256).reshape(1, 256)
cbar_ax.pcolormesh(bar_data, cmap=cmap)
fig.text(.29, .005 + cbar_height * .5, zthresh,
color="black", size=20, ha="right", va="center")
fig.text(.71, .005 + cbar_height * .5, round(zmax,2),
color="black", size=20, ha="left", va="center")
plt.savefig(filepath)
In [8]:
def calculate_sat_point(template, contrast, subj=None):
"""Calculate the point at which the colormap should saturate."""
data = []
for hemi in ["lh", "rh"]:
hemi_file = template.format(contrast=contrast, subj=subj, hemi=hemi)
hemi_data = nib.load(hemi_file).get_data()
data.append(hemi_data)
data = np.concatenate(data)
return max(3.71, moss.percentiles(data, 98))
In [9]:
def add_mask_overlay(b, mask_file):
"""Gray-out vertices outside of the common-space mask."""
mask_data = nib.load(mask_file).get_data()
# Plot the mask
mask_data = np.logical_not(mask_data.astype(bool)).squeeze()
if mask_data.any():
b.add_data(mask_data, min=0, max=10, thresh=.5,
colormap="bone", alpha=.6, colorbar=False)
In [35]:
def add_roi_overlay(b, roi_file, thresh, colormap, hemi):
"""Gray-out vertices outside of the common-space mask."""
reg_file = op.join(os.environ["FREESURFER_HOME"], "average/mni152.register.dat")
roi_vol = sf.project_volume_data(roi_file, hemi, reg_file)
# roi_vol = sf.project_volume_data(roi_file, hemi, subject_id='fsaverage')
# Plot the roi
if roi_vol.any():
b.add_data(roi_vol, min=0, max=0.05, thresh=thresh,
colormap=colormap, alpha=0.8, colorbar=False)
# b.add_contour_overlay(roi_vol, min=0, max=0.05, colormap = colormap,
# n_contours = 5, line_width=2.5, hemi=hemi)
# b.contour['colorbar'].visible = False
return roi_vol
In [11]:
def add_stat_overlay(b, stat_file, hemi, thresh, max, sign, sig_to_z=False):
"""Plot a surface-encoded statistical overlay."""
stat_data = nib.load(stat_file).get_data()
# Possibly convert -log10(p) images to z stats
if sig_to_z:
stat_sign = np.sign(stat_data)
p_data = 10 ** -np.abs(stat_data)
z_data = stats.norm.ppf(p_data)
z_data[np.sign(z_data) != stat_sign] *= -1
stat_data = z_data
# Plot the statistical data
plot_stat = ((sign == "pos" and (stat_data > thresh).any())
or (sign == "neg" and (stat_data < -thresh).any())
or (sign == "abs" and (np.abs(stat_data) > thresh).any()))
if plot_stat:
stat_data = stat_data.squeeze()
# b.add_overlay(stat_data, min=thresh, max=max, sign=sign, name="stat")
b.add_data(stat_data, min=thresh, max=max, colormap = 'OrRd_r',
thresh = 1e-5, alpha=.7, colorbar=False)
# b.add_contour_overlay(stat_data, min=thresh, max=max, colormap = 'OrRd_r',
# n_contours = 9, line_width=2, hemi=hemi)
return stat_data
In [12]:
def plot_contrast(subj, contrast, stat_temp, mask_temp, png_temp,
args, z_thresh, sign, viewpoints, hemis,
elevation=None, azimuths=None,
close_window=True, save_fig=True):
# Calculate where the overlay should saturate
z_max = calculate_sat_point(stat_temp, contrast, subj)
# save which should be manual
if 'manual' in viewpoints:
manual_index = viewpoints.index('manual')
# Create image panel
png_panes = []
for hemi in hemis:
# need a manual viewpoint for this hemi?
if azimuths is not None:
rot = dict(azimuth=azimuths[hemi],
elevation=elevation)
viewpoints[manual_index] = rot
print viewpoints
# Initialize the brain object
b_subj = subj if args.regspace == "epi" else "fsaverage"
b = Brain(b_subj, hemi, args.geometry,
config_opts=config_opts)
# Plot the mask
mask_file = mask_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
add_mask_overlay(b, mask_file)
# Plot the overlay
stat_file = stat_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
add_stat_overlay(b, stat_file, hemi, z_thresh, z_max, sign,
sig_to_z=args.regspace == "fsaverage")
# png_panes.append(view_panes(b, sign, hemi, viewpoints))
montage = b.save_montage(None, order = [viewpoints[0], 'ventral', 'medial'], orientation='v')[0:700,0:393,0:4]
print montage.shape
png_panes.append(montage)
# Close window if requested
if close_window:
b.close()
if save_fig:
# Stitch the hemisphere pngs together and save
full_png = np.hstack((png_panes[0], png_panes[1])) #np.concatenate(png_panes, axis=0) #1
png_file = png_temp.format(contrast=contrast, hemi=hemi, subj=subj)
print png_file
mplimg.imsave(png_file, full_png)
png_file = png_colorbar.format(contrast=contrast, hemi=hemi, subj=subj)
plot_colorbar(z_thresh, z_max, 'OrRd_r', png_file)
return png_panes
In [12]:
In [13]:
def view_panes(b, sign, hemi, viewpoints):
image_panes = []
for view in viewpoints:
# Handle the colorbar
if "stat" in b.overlays:
show_pos = True
show_neg = False
if sign in ("pos", "abs"):
b.overlays["stat"].pos_bar.visible = show_pos
if sign in ("neg", "abs"):
b.overlays["stat"].neg_bar.visible = show_neg
# Set the view and screenshot
# b.show_view(view, distance=330)
if view in ['med', 'ven']:
b.show_view(view, distance=330)
else:
b.show_view(view, distance=370)
sleep(0.2) # pause for show_view
image_panes.append(b.screenshot())
# Stitch the images together and return the array
full_image = np.concatenate(image_panes, axis=0)
return full_image
In [119]:
def plot_label(subj, label, mask_temp, alpha, color,
args, viewpoints, hemi, contrast, plot_contrast, b=None):
for view in viewpoints:
# Initialize the brain object
b_subj = subj if args.regspace == "epi" else "fsaverage"
if b is None:
b = Brain(b_subj, hemi, args.geometry,
config_opts=config_opts)
# Plot the mask
mask_file = mask_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
add_mask_overlay(b, mask_file)
if plot_contrast:
z_max = calculate_sat_point(stat_temp, contrast, subj)
stat_file = stat_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
add_stat_overlay(b, stat_file, hemi,
z_thresh, z_max, corr_sign,
sig_to_z=args.regspace == "fsaverage")
# Plot the label
b.add_label(label, alpha=alpha, color=color, borders=True)
b.show_view(view, distance=330)
return b
In [96]:
def plot_roi(subj, roi_file, mask_temp, alpha, thresh, colormap,
args, viewpoints, hemi, contrast, plot_contrast, b=None):
for view in viewpoints:
# Initialize the brain object
b_subj = subj if args.regspace == "epi" else "fsaverage"
if b is None:
b = Brain(b_subj, hemi, args.geometry,
config_opts=config_opts)
# Plot the mask
mask_file = mask_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
add_mask_overlay(b, mask_file)
# Plot the roi
roi_data = add_roi_overlay(b, roi_file, thresh, colormap, hemi)
if plot_contrast:
z_max = calculate_sat_point(stat_temp, contrast, subj)
stat_file = stat_temp.format(contrast=contrast,
hemi=hemi, subj=subj)
stat_data = add_stat_overlay(b, stat_file, hemi, z_thresh, z_max, corr_sign,
sig_to_z=args.regspace == "fsaverage")
else:
stat_data=[]
b.show_view(view, distance=330)
return b, stat_data, roi_data
In [90]:
contrasts
Out[90]:
In [ ]:
mlab.view()
In [26]:
contrast = 'hit-miss'
elevation = 82
azimuth = {"lh": -134}
hemi = ["lh"]
# 134 lh, 46 rh (90 degrees +- 44)
viewpoints = ['manual']
plot_contrast(subj, contrast, stat_temp, mask_temp, png_temp,
args, z_thresh, corr_sign, viewpoints, hemi,
elevation=elevation, azimuths=azimuth,
close_window=False, save_fig=False)
In [94]:
contrast = 'cr-2_fa'
viewpoints = ['manual', 'ven', 'med']
# viewpoints = ['lat', 'ven', 'med']
azimuths = None
elevation = None
# specify if a viewpoint is 'manual'
elevation = 82
azimuths = {'lh': -134, 'rh': -46}
hemis = ["lh", "rh"]
png_panes = plot_contrast(subj, contrast, stat_temp, mask_temp, png_temp,
args, z_thresh, corr_sign,
viewpoints, hemis, elevation=elevation, azimuths=azimuths)
In [95]:
plt.imshow(np.hstack((png_panes[0], png_panes[1])))
Out[95]:
In [129]:
label = 'lateraloccipital'
save_file = True
color= 'navy'
alpha = 1
# rot = dict(azimuth=-134, elevation=82)
rot = dict(azimuth=-46, elevation=82)
# rot = 'med'
viewpoints = [rot]
hemi = 'rh'
contrast = 'study_rep'
b = plot_label(subj, label, mask_temp, alpha, color,
args, viewpoints, hemi, contrast, plot_contrast=True)
if save_file:
png_file = '/Users/steph-backup/Dropbox/Experiments/ObjFam/study_rep_lo_rh.png'
b.save_image(png_file)
In [96]:
%matplotlib inline
import seaborn as sns
In [119]:
Out[119]:
In [104]:
sns.palplot(sns.color_palette('OrRd_r', n_colors = 50))
In [68]:
import seaborn.apionly as sns
sns.choose_colorbrewer_palette("s");
In [71]:
x = np.random.rand(30, 30)
In [76]:
sns.heatmap(x, cmap=cmap)
Out[76]:
In [39]:
label = '17Networks_5'; color = 'hotpink'
alpha = .5
# viewpoints = ['lat']
rot = dict(azimuth=-134, elevation=82)
viewpoints = [rot]
hemi = 'lh'
contrast = 'study_rep'
b = plot_label(subj, label, mask_temp, alpha, color,
args, viewpoints, hemi, contrast,
plot_contrast=False)
label = 'Network5_SPL'
color = 'orange'
b = plot_label(subj, label, mask_temp, alpha, color,
args, viewpoints, hemi, contrast,
plot_contrast=False)
save_file = False
if save_file:
png_file = '/Users/steph-backup/Experiments/ObjFam/data/fsaverage/17Networks_5.png'
b.save_image(png_file)
In [67]:
png_file = png_label_temp.format(contrast=contrast, hemi=hemi, subj=subj, label=label)
b.save_image(png_file)
In [111]:
roi_file = '/Users/steph-backup/Experiments/ObjFam/data/fsaverage/masks/topdown_ALE_pN05_fslMNI.nii.gz'
roi = 'topdown_ALE_pN05'
alpha = 0.4
rot = dict(azimuth=-134, elevation=82)
viewpoints=[rot]
hemi = 'lh'
contrast = '1_old'
# contrast = 'linear-study'
# colormap = 'PuBuGn_r'
colormap = 'Greens_r'
# colormap=['darkgreen', 'darkgreen']
b, stat_data, roi_data = plot_roi(subj, roi_file, mask_temp, alpha, 0.00001, colormap,
args, viewpoints, hemi, contrast,
plot_contrast=False)
# roi_file = '/Users/steph-backup/Experiments/ObjFam/data/fsaverage/masks/topdownALEventral_fslMNI_thresh.nii.gz'
# colormap=['lawngreen', 'lawngreen']
# b, stat_data, roi_data = plot_roi(subj, roi_file, mask_temp, alpha, 0.001, colormap,
# args, viewpoints, hemi, contrast,
# plot_contrast=False, b=b)
coords = [[-25.1, -61.07,51.19],
[-29.32,-58.13,59.93],
[-31.6,-48.12,52.26],
[-16.39,-62.75,55.69],
[-26.12,-54.07,57.24],
[-33.85,-42.17,47.23],
[-27.57,-78.67,23.84],
[-25.16,-52.77,48.13],
[-27.46,-75.71,32.51],
[-25.03,-64.91,56.05],
[-25.1,-54.48,52.78],
[-23.09,-66.57,38.25],
[-25.05,-51.94,57.01],
[-18.4,-59.39,68.84],
[-25.36,-76.12,28.03],
[-29.33,-59.3,58.93],
[-31.82,-80.4,28.56],
[-17.16,-65.89,79.55],
[-44.65,-35.62,49],
[-23,-79,32],
[-12,-44,28],
[-44,-64,32],
[-44,-48,36],
[-16,-52,56],
[-16,-56,56],
[-40,-64,28],
[-40,-48,32],
[-30.48,-55.46,54.09],
[-37,-32,59],
[-27,-78,38],
[-18,-63,54],
[-39.12,-46.67,56.73],
[-22.83,-53.64,61.62],
[-12.06,-56.05,58.31],
[-39.17,-47.09,52.29],
[-27.42,-57.86,39.71],
[-18.52,-57.03,59.64],
[-21.88,-59.04,49.81],
[-39.3,-40.37,43.78],
[-22,-62,60],
[-42.35,-53.17,56.3],
[-20.75,-75.09,50.26],
[-17.39,-68.72,60.78],
[-35.94,-52.6,50.54],
[-25.37,-73.99,27.82],
[-29.28,-53.46,63.95],
[-31.57,-50.04,54.69],
[-29.54,-46.61,45.35],
[-11.95,-64.03,64.71],
[-22.88,-69.37,54.21],
[-12.05,-75.71,54.66],
[-29.62,-84.44,31.17],
[-29.28,-53.46,63.95],
[-29.54,-44.49,45.14],
[-25.37,-71.87,27.61],
[-22.86,-67.03,56.22],
[-16.51,-70.17,45.22],
[-22.96,-75.41,46.97],
[-18.51,-73.4,56.78],
[-23,-73,40],
[-21,-68,52],
[-25,-61,55],
[-26,-57,54],
[-33.72,-53.24,55.04],
[-19.63,-68.19,55.16],
[-30.66,-31.62,45.01],
[-22,-72,34],
[-18.62,-61.05,51.08],
[-27.24,-52.16,54.83],
[-25.31,-73.58,32.26],
[-31.6,-49.19,52.37],
[-29.65,-79.22,29.53],
[-26.32,-62.11,40.11],
[-14.24,-61.68,55.55]]
b.add_foci(coords, map_surface="white", color="white", scale_factor=.5, alpha=1)
In [112]:
png_file = png_roi_temp.format(contrast=contrast, hemi=hemi, subj=subj, roi=roi)
print png_file
png_file = '/Users/steph-backup/Dropbox/Experiments/ObjFam/alemap_foci.png'
b.save_image(png_file)
In [ ]:
label_atten = '17Networks_5'; color = 'dodgerblue'
label_medialpar = '17Networks_11'; color = 'mediumseagreen'
label_frontoparietal = '17Networks_12'; color = 'purple'
label_default = '17Networks_16'; color = 'cornflowerblue'
In [ ]:
In [81]:
brain = Brain('subj03', 'lh', "inflated", views=dict(azimuth=-134, elevation=82),
config_opts=dict(background="white"))
mri_file = "/Users/steph-backup/Experiments/ObjFam/data/subj03/masks/lh.17Networks_5_SPL.nii.gz"
reg_file = "/Users/steph-backup/Experiments/ObjFam/scripts/analysis/objfam/subj03/preproc/run_1/func2anat_tkreg.dat"
surf_data_hemi = io.project_volume_data(mri_file, "lh", reg_file=reg_file)
brain.add_data(surf_data_hemi, hemi='lh', colorbar=False, min=0, max=1, colormap=['darkorchid', 'darkorchid'], thresh=.1, alpha=.7)
mri_file = "/Users/steph-backup/Experiments/ObjFam/data/subj03/masks/lh.17Networks_5_LO.nii.gz"
reg_file = "/Users/steph-backup/Experiments/ObjFam/scripts/analysis/objfam/subj03/preproc/run_1/func2anat_tkreg.dat"
surf_data_hemi = io.project_volume_data(mri_file, "lh", reg_file=reg_file)
brain.add_data(surf_data_hemi, hemi='lh', colorbar=False, min=0, max=1, colormap=['deeppink', 'deeppink'], thresh=.1, alpha=.7)
png_file = '/Users/steph-backup/Dropbox/Experiments/ObjFam/17Networks_5_SPL_LO.png'
brain.save_image(png_file)
In [ ]: